#'
#' Calculate various topographic variables based on elevation
#'
#'
#' @param topo a dataframe contains three columns named "x", "y" and "elev"
#' @param plotdim the size of the whole plot
#' @param gridsize_org the original scale of the given topo variable
#' @param edgecorrect used in adjust the convexity in the edge
#' 
#' 
#' @details 
#' it will calculate the mean elevation, convexity, slope, aspect, 
#' wetness index and altitude above channel network at the scale of gridesize*4
#' 
#' @note 
#' it needs adehabitat, rgdal and RSAGA packages. it also need the SAGA-GIS software 
#' (http://sourceforge.net/projects/saga-gis/files/SAGA%20-%202.1/) installed under R package 
#' file "RSAGA" in R library.
#'

#' @export
topo_variables=function(topo,plotdim,gridsize_org=5,edgecorrect=TRUE){
  if(!inherits(topo,"env"))
    stop("topo should be env object")
  #define the gride size in topo object
  topo=as.data.frame(topo)
  names(topo)=c("x","y","elev")
  gridsize=gridsize_org*4
  elev.list=elev.to.list(topo,gridsize=gridsize_org)
  topovb=allquadratslopes(elev.list,gridsize=gridsize,plotdim=plotdim,edgecorrect=edgecorrect)
  aspect=allquadrataspect(elev.list,gridsize=gridsize,plotdim=plotdim)
  topovb$aspect=aspect[,3]
  
  require(adehabitat)
  require(rgdal)
  require(RSAGA)
  
  elev=subset(topovb,select="meanelev")
  twi_ach=twi_ach(elevation=elev,plotdim=plotdim,gridsize=gridsize)
  topovb=cbind(topovb,twi_ach)
  return(topovb)
}



elev.to.list=function(elevfile,gridsize)
{ 
  xdim=max(elevfile$x)
  ydim=max(elevfile$y)
  elevmat=matrix(elevfile$elev,nrow=1+ydim/gridsize,ncol=1+xdim/gridsize,byrow=F)
  return(list(col=elevfile,mat=elevmat))
}

quadslope=function(cornerelev,gridsize)
{
  slope=numeric(4)
  z=numeric(3)
  for(j in 1:4)
  {
    post=1
    for(k in (j+1):(j+3))
    {
      if(k>4) m=k%%4
      else m=k
      z[post]=cornerelev[m]
      post=post+1
    }
    slope[j]=calcslope(z,gridsize)
  }
  return( c(mean(slope),sqrt(var(slope))) )
}

calcslope=function(z, gridsize)
{
  z2=z[3]-z[2]
  z1=z[1]-z[2]
  if(z2==0)
  {
    if(z1==0) return(0)
    else denom=sqrt( 1+(gridsize/z1)^2 )
    theta1 = acos( (-gridsize/z1)/denom )
    theta2 = acos( (gridsize/z1)/denom )
  }
  else 
  {
    denom = sqrt( 1+(z1/z2)^2+(gridsize/z2)^2 )
    theta1 = acos( (-gridsize/z2)/denom )
    theta2 = acos( (gridsize/z2)/denom )
  }
  if(theta1<=theta2) return(180*theta1/pi)
  else return(180*theta2/pi)
}

rowcol.to.index=function(rowno,colno,gridsize,plotdim)
{
  badrc=(rowno<=0 | colno<=0 | rowno>plotdim[2]/gridsize | colno>plotdim[1]/gridsize)
  rowno=rowno-1
  colno=colno-1
  maxrow=floor(plotdim[2]/gridsize)
  index=colno*maxrow+rowno+1
  if(length(badrc[badrc>0])) index[badrc]=NA
  return(index)
}

findborderquads=function(index,dist,gridsize,plotdim)
{
  bound.index=numeric(8)
  no.boundaries=0
  row=index.to.rowcol(index,gridsize,plotdim)$row
  col=index.to.rowcol(index,gridsize,plotdim)$col
  maxrow=plotdim[2]/gridsize
  maxcol=plotdim[1]/gridsize
  layers=floor(dist/gridsize)
  for(i in (row-layers):(row+layers))
    for(j in (col-layers):(col+layers))
      if(i!=row | j!=col)
        if(i>=1 & i<=maxrow & j>=1 & j<=maxcol)
        {
          no.boundaries=no.boundaries+1
          bound.index[no.boundaries]=rowcol.to.index(i,j,gridsize,plotdim)
        }
  return( bound.index[bound.index>0] )
} 

index.to.gxgy=function(index,gridsize,plotdim)
{
  badindex=(index<=0 | index>plotdim[1]*plotdim[2]/(gridsize^2))
  rc=index.to.rowcol(index,gridsize,plotdim)
  gx=gridsize*(rc$col-1)
  gy=gridsize*(rc$row-1)
  if(length(badindex[badindex>0])) gx[badindex]=gy[badindex]=-1
  return(data.frame(gx=gx,gy=gy))
}

index.to.rowcol=function(index,gridsize,plotdim)
{
  index=index-1
  badindex=(index<0 | index>=plotdim[1]*plotdim[2]/(gridsize^2))
  maxrow=floor(plotdim[2]/gridsize)
  rowno=index%%maxrow
  colno=floor((index-rowno)/maxrow)
  row=rowno+1
  col=colno+1
  if(length(badindex[badindex>0])) row[badindex]=col[badindex]=-1
  return(data.frame(row=row,col=col))
}

allquadratslopes=function(elev,gridsize,plotdim, edgecorrect=TRUE)
{
  rw=cl=0
  on.exit(cat(rw," ",cl,"\n"))
  columns=1+max(elev$col$x)/gridsize
  rows=1+max(elev$col$y)/gridsize
  totalquads=(columns-1)*(rows-1)
  cat("Calculating topographic indices for ", totalquads, " quadrats\n")
  elevdata=elev$col[elev$col$x%%gridsize==0 & elev$col$y%%gridsize==0,]
  elevmat=matrix(elevdata$elev,nrow=rows,ncol=columns,byrow=F)
  meanelev=convex=convex2=slope=numeric()
  corner=sideht=numeric()
  # Mean elevation of four corners
  for(c in 1:(columns-1))
    for(r in 1:(rows-1))
    {
      quad.index=rowcol.to.index(r,c,gridsize=gridsize,plotdim=plotdim)
      corner[1]=elevmat[r,c]
      corner[2]=elevmat[r+1,c]
      corner[3]=elevmat[r+1,c+1]
      corner[4]=elevmat[r,c+1]
      meanelev[quad.index]=mean(corner)
      slope[quad.index]=quadslope(corner,gridsize=gridsize)[1]
      if(c%%33==0 & r%%33==0) cat("Finding elevation and slope of quadrat ", quad.index, "\n")
    }
  # Convexity
  for(i in 1:totalquads)
  { 
    neighbor.quads=findborderquads(i,dist=gridsize,gridsize=gridsize,plotdim=plotdim)
    meanelev.neighbor=mean(meanelev[neighbor.quads])
    convex[i]=meanelev[i]-meanelev.neighbor
    if(i%%1000==0) cat("Finding convexity of quadrat ", i, "\n")
  }
  # correcting convexity in edge quadrats, based on center of the 20x20 rather
  # than surrounding 20x20s. This requires that the elev$mat has an elevation
  # at the middle of every grid cell.
  if(edgecorrect)
  {
    for(c in 1:(columns-1))
      for(r in 1:(rows-1))
      {
        if((c==1) | (c==(columns-1)) | (r==1) | (r==(rows-1)))
        { 
          quad.index=rowcol.to.index(r,c,gridsize=gridsize,plotdim=plotdim)
          xy=index.to.gxgy(quad.index,gridsize=gridsize,plotdim=plotdim)
          midx=xy$gx+gridsize/2
          midy=xy$gy+gridsize/2
          #browser()
          midelev=subset(elev$col,x==midx & y==midy)$elev
          convex[quad.index]=midelev-meanelev[quad.index]
        }
      }
  }
  return(data.frame(meanelev=meanelev,convex=convex,slope=slope))
}	

#calculating aspect for each quadrat
allquadrataspect=function(elev,gridsize,plotdim)
{
  rw=cl=0
  on.exit(cat(rw," ",cl))
  columns=1+max(elev$col$x)/gridsize
  rows=1+max(elev$col$y)/gridsize
  totalquads=(columns-1)*(rows-1)
  elevdata=elev$col[elev$col$x%%gridsize==0 & elev$col$y%%gridsize==0,]
  elevmat=matrix(elevdata$elev,nrow=rows,ncol=columns,byrow=F)
  xmat=matrix(elevdata$x,nrow=rows,ncol=columns,byrow=F)
  ymat=matrix(elevdata$y,nrow=rows,ncol=columns,byrow=F)
  aspect=x.co=y.co=numeric()
  corner=xcorner=ycorner=numeric()
  for(c in 1:(columns-1)) 
    for(r in 1:(rows-1))
    {
      quad.index=rowcol.to.index(r,c,gridsize=gridsize,plotdim=plotdim)
      corner[1]=elevmat[r,c];corner[2]=elevmat[r+1,c];corner[3]=elevmat[r+1,c+1];corner[4]=elevmat[r,c+1]
      xcorner[1]=xmat[r,c]; xcorner[2]=xmat[r+1,c]; xcorner[3]=xmat[r+1,c+1]; xcorner[4]=xmat[r,c+1]
      ycorner[1]=ymat[r,c]; ycorner[2]=ymat[r+1,c]; ycorner[3]=ymat[r+1,c+1]; ycorner[4]=ymat[r,c+1]
      x.co[quad.index]=mean(xcorner);y.co[quad.index]=mean(ycorner)   
      aspect[quad.index]=quadaspect(corner,gridsize=gridsize)  
    }
  return(data.frame(x.co,y.co,aspect=aspect))
}

quadaspect=function(cornerelev,gridsize)
{
  aspect=NA
  fy=cornerelev[2]+cornerelev[3]-cornerelev[1]-cornerelev[4]
  fx=cornerelev[3]+cornerelev[4]-cornerelev[1]-cornerelev[2]
  if(fx==0){
    if(fy<0) aspect=0
    else if(fy>0) aspect=180
  }
  else aspect=180-atan(fy/fx)*180/pi+90*fx/abs(fx)
  if(is.na(aspect)) browser();
  return(aspect)
}

#how to get SAGA-GIS 2.1.4 software: http://sourceforge.net/projects/saga-gis/files/SAGA%20-%202.1/, and you should install it
#under R package file "RSAGA" in R library(e.g. "C:/softwares/R-3.1.2/library/RSAGA/SAGA-GIS").
#elevation: elevation is one-column data, meanelev, which is the mean elevation of all quadrats
#plotdim: is the dimension of the entire plot(e.g. plotdim(bci)=c(1000,500))
#gridsize: is the size of each quadrat
twi_ach=function(elevation,plotdim,gridsize)
{
  #transform elevation data to ASCII raster file
  elev.asc=as.asc(t(matrix(elevation$meanelev,plotdim[2]/gridsize,plotdim[1]/gridsize)),xll=gridsize/2,yll=gridsize/2,gridsize)
  export.asc(elev.asc,"elev.asc")
  
  #converts grid files from ESRI's ASCII (.asc) format to SAGA's grid format (.sgrd).
  rsaga.esri.to.sgrd(in.grids="elev.asc", out.sgrds="elev.sgrd", in.path=getwd())
  
  #calculate topographic wetness index using SAGA-GIS
  rsaga.wetness.index(in.dem="elev.sgrd",out.wetness.index="twi.sgrd")
  
  # catchment area:
  rsaga.geoprocessor(lib="ta_preprocessor", module=2, param=list(DEM="elev.sgrd", DEM_PREPROC="sinks.sgrd"))
  # extract channel network:
  rsaga.geoprocessor(lib="ta_channels", module=0, param=list(ELEVATION="elev.sgrd", CHNLNTWRK="chnw.sgrd",INIT_GRID="sinks.sgrd"))
  # altitude above channel network:
  rsaga.geoprocessor(lib="ta_channels", module=3, param=list(ELEVATION="elev.sgrd", CHANNELS="chnw.sgrd", DISTANCE="achan.sgrd", THRESHOLD=0.1, NOUNDERGROUND=TRUE))#
  
  # converts grid files from SAGA's grid format (.sgrd) format to ESRI's ASCII (.asc). 
  rsaga.sgrd.to.esri(in.sgrds=c("twi.sgrd", "achan.sgrd"), out.grids=c("twi.asc", "achan.asc"), prec=2, out.path=getwd())
  #get topographic wetness index
  twi_matrix = read.ascii.grid("twi.asc", return.header=F)
  twi_matrix2 = cbind(as.matrix(1:(plotdim[2]/gridsize), plotdim[2]/gridsize,1), twi_matrix)
  twi = as.numeric(twi_matrix2[order(twi_matrix2[,1], decreasing = T),][,-1])
  #get vertical distance from the channel network 
  achan_matrix = read.ascii.grid("achan.asc", return.header = F)
  achan_matrix2 = cbind(as.matrix(1:(plotdim[2]/gridsize), plotdim[2]/gridsize, 1), achan_matrix)
  achan=as.numeric(achan_matrix2[order(achan_matrix2[,1], decreasing = T),][,-1])
  return(data.frame(twi=twi,achan=achan))
}	
